#!/usr/bin/perl

# ------------------------------------------------------------------------------
# --- Get homologous sequences from CD-HIT output so we can duplicate final
# ---  MySQL output
# ------------------------------------------------------------------------------

use strict;
use warnings;

my $usage = "Usage: perl $0\n";

my $in_clstr = "pdb_dssp_seq_95.fa.clstr";

my $in_sql = "combined.sql";
my $out_sql = $in_sql . ".homolog.sql";

# Percent ID cut-off. To not include homologs with < 95% similarity, for example,
# set $cutoff = 95
my $cutoff = 100;

open (CLST, "<$in_clstr")
	or die "ERROR: Could not open input cluster file, $in_clstr\n";
open (OUT, ">$out_sql")
	or die "ERROR: Could not open output sql file, $out_sql\n";

# Run though cluster file and figure out the starred PDB and chain for each cluster

my $clust_num;

my %starred_pdb;
my %starred_chain;

while (my $line = <CLST>) {

	if ($line =~ /^>Cluster (\d+)/) {
		$clust_num = $1;
	} elsif ($line =~ /\*/) {
		if ($line =~ /\d+\t\d+aa, >(.{4}):(.):sequence.../) {
			$starred_pdb{$clust_num} = $1;
			$starred_chain{$clust_num} = $2;
		}
	}
}

close CLST;

# Run though cluster file again, but this time grab the relevant SQL statements and 
# swap in the homologous PDB and chain

open (CLST, "<$in_clstr")
	or die "ERROR: Could not open input cluster file, $in_clstr\n";

my @grep_result_1;
my @grep_result_2;

my $hom_pdb;
my $hom_chain;

while (my $line = <CLST>) {

	if ($line =~ /^>Cluster (\d+)/) {
		$clust_num = $1;
	
		print STDERR "Parsing Cluster #" . $clust_num . "...\n";
		
		$hom_pdb = $starred_pdb{$clust_num};
		$hom_chain = $starred_chain{$clust_num};
		
		# Find lines that match this PDB and chain in SQL output
		# ...with matching chain in first position:
		my $grep_cmd_1 = 
			'grep "VALUES (\"' . $hom_pdb . '\", \"' . $hom_chain . '" '   . $in_sql;
		# ...with matching chain in second position:
		my $grep_cmd_2 = 
			'grep "VALUES (\"' . $hom_pdb . '\", \". ' . $hom_chain . '" ' . $in_sql;
		
		@grep_result_1 = ();
		@grep_result_2 = ();
				
		@grep_result_1 = `$grep_cmd_1`;
		@grep_result_2 = `$grep_cmd_2`;
		
	} else {
		if ($line =~ /\d+\t\d+aa, >(.{4}):(.):sequence... at ([\d\.]+)%/) {
			my $pair_pdb = $1;
			my $pair_chain = $2;
			
			my $perc_match = $3;
			if ($perc_match < $cutoff) {
				print "\tSkipping $pair_pdb $pair_chain. % Match = [$perc_match]...\n";
			} else {
				# Replace original with homolog
				# First those matching chain in first position:
				foreach my $result_1 (@grep_result_1) {
					my $new_sql = $result_1;
					$new_sql =~ s/(VALUES \(\")($hom_pdb)/$1$pair_pdb/g;
					$new_sql =~ s/(VALUES \(\".{4}\", \")./$1$pair_chain/g;
					print OUT $new_sql;
				}

				# Then those matching chain in second position:
				foreach my $result_2 (@grep_result_2) {
					my $new_sql = $result_2;
					$new_sql =~ s/(VALUES \(\")($hom_pdb)/$1$pair_pdb/g;
					$new_sql =~ s/(VALUES \(\".{4}\", \". )./$1$pair_chain/g;
					print OUT $new_sql;
				}	
			}
		} 
	}
}

exit;